library(dplyr)
library(stringr)
library(multiwayvcov)
library(lmtest)
library(ggplot2)
library(cregg)
library(stargazer)
library(margins)

setwd("C:/Users/jhpar/OneDrive/Desktop/Research/P4P/Data/SecondPaper")
load(file='P4P_Clean.rda')
str(P4P_Clean)

P4P_Clean <- within(P4P_Clean, `Total pay` <- relevel(`Total pay`, ref = "Slightly below average"))
P4P_Clean <- within(P4P_Clean, `Performance bonuses` <- relevel(`Performance bonuses`, ref = "No performance bonuses; fixed salary"))
P4P_Clean <- within(P4P_Clean, `Performance bonuses binary` <- relevel(`Performance bonuses binary`, ref = "No (fixed salary)"))
P4P_Clean <- within(P4P_Clean, `Job performance evaluation` <- relevel(`Job performance evaluation`, ref = "A supervisor evaluation of your work"))
P4P_Clean <- within(P4P_Clean, `Goal based evaluation` <- relevel(`Goal based evaluation`, ref = "No (supervisor based)"))
P4P_Clean <- within(P4P_Clean, `Current community involvement` <- relevel(`Current community involvement`, ref = "Rare participation"))
P4P_Clean <- within(P4P_Clean, `Community income` <- relevel(`Community income`, ref = "Low income"))
P4P_Clean <- within(P4P_Clean, `Community demographics` <- relevel(`Community demographics`, ref = "Mostly white"))
P4P_Clean <- within(P4P_Clean, `Overtime work` <- relevel(`Overtime work`, ref = "Never required"))
P4P_Clean <- within(P4P_Clean, `Key job task` <- relevel(`Key job task`, ref = "Analysis identifying community needs"))
P4P_Clean <- within(P4P_Clean, PSM_C <- relevel(PSM_C, ref = "Low"))
P4P_Clean <- within(P4P_Clean, Efficacy_C <- relevel(Efficacy_C, ref = "Low"))

P4P_Clean_Binary <- P4P_Clean %>%
  mutate(`Total pay binary`=as.factor(ifelse(`Total pay`=="Slightly below average","Slightly below average","About average and Slightly above average"))) %>%
  mutate(`Current community involvement binary`=as.factor(ifelse(`Current community involvement`=="Rare participation","Rare participation","Moderate and Frequent participation"))) %>%
  mutate(`Community income binary`=as.factor(ifelse(`Community income`=="Low income","Low income","Average and High income"))) %>%
  mutate(`Community demographics binary`=as.factor(ifelse(`Community demographics`=="Mostly white","Mostly white","Non-White"))) %>%
  mutate(`Overtime work binary`=as.factor(ifelse(`Overtime work`=="Never required","Never required","Occasionally and Frequently required"))) %>%
  mutate(`Key job task binary`=as.factor(ifelse(`Key job task`=="Analysis identifying community needs","Analysis identifying community needs","Interaction with community and organizational members")))
P4P_Clean_Binary <- within(P4P_Clean_Binary, `Total pay binary` <- relevel(`Total pay binary`, ref = "Slightly below average"))
P4P_Clean_Binary <- within(P4P_Clean_Binary, `Current community involvement binary` <- relevel(`Current community involvement binary`, ref = "Rare participation"))
P4P_Clean_Binary <- within(P4P_Clean_Binary, `Community income binary` <- relevel(`Community income binary`, ref = "Low income"))
P4P_Clean_Binary <- within(P4P_Clean_Binary, `Community demographics binary` <- relevel(`Community demographics binary`, ref = "Mostly white"))
P4P_Clean_Binary <- within(P4P_Clean_Binary, `Overtime work binary` <- relevel(`Overtime work binary`, ref = "Never required"))
P4P_Clean_Binary <- within(P4P_Clean_Binary, `Key job task binary` <- relevel(`Key job task binary`, ref = "Analysis identifying community needs"))


# 1. Baseline Model
# 1.1 Marginal Means
# Note: Since the 'cregg' package does not allow for two-way clustering,
# we cannot use the 'cregg'. 
# While the 'cregg' uses glm() and svgglm() to estimate marginal means,
# we use lm() to do so. Differences in standard errors are negligible. 
Base_tp <- lm(Chosen_Job ~ 0 + `Total pay binary`, data=P4P_Clean_Binary)
Base_tp.cl <- cluster.vcov(Base_tp, cbind(P4P_Clean_Binary$ID,P4P_Clean_Binary$Treatment.ID))
tp <- broom::tidy(coeftest(Base_tp, Base_tp.cl))

Base_pb <- lm(Chosen_Job ~ 0 + `Performance bonuses binary`, data=P4P_Clean_Binary)
Base_pb.cl <- cluster.vcov(Base_pb, cbind(P4P_Clean_Binary$ID,P4P_Clean_Binary$Treatment.ID))
pb <- broom::tidy(coeftest(Base_pb, Base_pb.cl))

Base_ge <- lm(Chosen_Job ~ 0 + `Goal based evaluation`, data=P4P_Clean_Binary)
Base_ge.cl <- cluster.vcov(Base_ge, cbind(P4P_Clean_Binary$ID,P4P_Clean_Binary$Treatment.ID))
ge <- broom::tidy(coeftest(Base_ge, Base_ge.cl))

Base_cci <- lm(Chosen_Job ~ 0 + `Current community involvement binary`, data=P4P_Clean_Binary)
Base_cci.cl <- cluster.vcov(Base_cci, cbind(P4P_Clean_Binary$ID,P4P_Clean_Binary$Treatment.ID))
cci <- broom::tidy(coeftest(Base_cci, Base_cci.cl))

Base_ci <- lm(Chosen_Job ~ 0 + `Community income binary`, data=P4P_Clean_Binary)
Base_ci.cl <- cluster.vcov(Base_ci, cbind(P4P_Clean_Binary$ID,P4P_Clean_Binary$Treatment.ID))
ci <- broom::tidy(coeftest(Base_ci, Base_ci.cl))

Base_cd <- lm(Chosen_Job ~ 0 + `Community demographics binary`, data=P4P_Clean_Binary)
Base_cd.cl <- cluster.vcov(Base_cd, cbind(P4P_Clean_Binary$ID,P4P_Clean_Binary$Treatment.ID))
cd <- broom::tidy(coeftest(Base_cd, Base_cd.cl))

Base_ow <- lm(Chosen_Job ~ 0 + `Overtime work binary`, data=P4P_Clean_Binary)
Base_ow.cl <- cluster.vcov(Base_ow, cbind(P4P_Clean_Binary$ID,P4P_Clean_Binary$Treatment.ID))
ow <- broom::tidy(coeftest(Base_ow, Base_ow.cl))

Base_kt <- lm(Chosen_Job ~ 0 + `Key job task binary`, data=P4P_Clean_Binary)
Base_kt.cl <- cluster.vcov(Base_kt, cbind(P4P_Clean_Binary$ID,P4P_Clean_Binary$Treatment.ID))
kt <- broom::tidy(coeftest(Base_kt, Base_kt.cl))

Base_MM <- rbind(tp,pb,ge,cci,ci,cd,ow,kt) %>%
  mutate(upper=estimate+(qnorm(0.975)*std.error), lower=estimate-(qnorm(0.975)*std.error)) 


# 1.2 Plot (Figure 2)
Base.Cat <- data.frame(
  Value=c("Total pay:", "Slightly below average","About average and Slightly above average",
          "Performance bonuses binary:", "No (fixed salary)","Yes (bonuses)",
          "Goal based evaluation:","No (supervisor based)","Yes (goal based)",
          "Current community involvement:","Rare participation","Moderate and Frequent participation",
          "Community income:","Low income","Average and High income",
          "Community demographics:","Mostly white","Non-White",
          "Overtime work:","Never required","Occasionally and Frequently required",
          "Key job task:","Analysis identifying community needs","Interaction with community and organizational members"))

Base_MM.Match <- data.frame(trimws(gsub(".*[`]","",Base_MM$term),"both"))
names(Base_MM.Match) <- "Value"
Base_MM.Match1 <- cbind(Base_MM,Base_MM.Match)

Base.MM.table <- merge(x=Base.Cat, y=Base_MM.Match1,
                       by='Value', all.x=TRUE)

Base.MM.table$Value <- factor(Base.MM.table$Value,
                              levels=c("Total pay:", "Slightly below average","About average and Slightly above average",
                                       "Performance bonuses binary:", "No (fixed salary)","Yes (bonuses)",
                                       "Goal based evaluation:","No (supervisor based)","Yes (goal based)",
                                       "Current community involvement:","Rare participation","Moderate and Frequent participation",
                                       "Community income:","Low income","Average and High income",
                                       "Community demographics:","Mostly white","Non-White",
                                       "Overtime work:","Never required","Occasionally and Frequently required",
                                       "Key job task:","Analysis identifying community needs","Interaction with community and organizational members"))

ggplot(Base.MM.table, aes(x = estimate, y = Value)) + 
  geom_vline(aes(xintercept = 0.5), size = .25, linetype = "dashed") + 
  geom_errorbarh(data = Base.MM.table, aes(xmax = upper, xmin = lower), size = .3, height = .6, position = position_nudge(y = 0), color="black") +
  geom_point(data = Base.MM.table, shape=19, size = 1.5, color='black') +
  scale_y_discrete(limits=rev) +
  theme(axis.text.y = element_text(hjust=0)) +
  theme(axis.text.y = element_text(face = 
                                     c("plain","plain","bold","plain","plain","bold","plain","plain",
                                       "bold","plain","plain","bold","plain","plain","bold","plain","plain","bold","plain","plain","bold",
                                       "plain","plain","bold"))) +
  ylab("") +
  xlab("Marginal Mean") +
  theme(axis.text.y = element_text(size=12, colour = "black")) +
  theme(text = element_text(size=12, colour = "black"))

## Comparable point estimates: One-way clustering!
mm(P4P_Clean_Binary,
   Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
     `Current community involvement binary` + `Community income binary` + 
     `Community demographics binary` + `Overtime work binary` + `Key job task binary`,
   id = ~ID)


# 1.3 AMCE (Appendix Table D1)
Base.AMCE <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                  `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                  `Overtime work binary` + `Key job task binary`, 
                data=P4P_Clean_Binary)
summary(Base.AMCE)
Base.AMCE.cl <- cluster.vcov(Base.AMCE, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
Base.stat <- coeftest(Base.AMCE, Base.AMCE.cl)
Twoway.se.Base <- sqrt(diag(Base.AMCE.cl))

# BH Procedure
pval <- Base.stat[,4]
p.bh <- p.adjust(pval, "fdr")
p.bh05 <- ifelse(p.bh<=0.05, 1, 0)
p.bh10 <- ifelse(p.bh<=0.10, 1, 0)


Base.AMCE.Collapse <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Job performance evaluation` +
                           `Current community involvement` + `Community income` + `Community demographics` +
                           `Overtime work` + `Key job task`, 
                         data=P4P_Clean_Binary)
summary(Base.AMCE.Collapse)
Base.AMCE.Collapse.cl <- cluster.vcov(Base.AMCE.Collapse, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
Base.Collapse.stat <- coeftest(Base.AMCE.Collapse, Base.AMCE.Collapse.cl)
Twoway.se.Base.Collapse <- sqrt(diag(Base.AMCE.Collapse.cl))

# BH Procedure
pval <- Base.Collapse.stat[,4]
p.bh <- p.adjust(pval, "fdr")
p.bh05 <- ifelse(p.bh<=0.05, 1, 0)
p.bh10 <- ifelse(p.bh<=0.10, 1, 0)


stargazer(Base.AMCE, Base.AMCE.Collapse,
          se=list(Twoway.se.Base,Twoway.se.Base.Collapse), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          order = c(1,9,10,2,11,12,13,3,14,15,16,4,17,18,5,19,20,6,21,22,23,7,24,25,8,26,27,28),
          covariate.labels=c("Total pay: About average and Slightly above average",
                             "Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%)",
                             "Performance bonuses: moderate (10%)",
                             "Performance bonuses: small (5%)",
                             "Goal based evaluation: Yes (B)",
                             "Goal based evaluation: Attendance",
                             "Goal based evaluation: Changes",
                             "Goal based evaluation: Satisfaction",
                             "Current community involvement binary: Moderate and Frequent participation",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income binary: Average and High income",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics binary: Non-White",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work binary: Occasionally and Frequently required",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task binary: Interaction with community and organizational members",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors"),
          out="Appendix_Table_D1.txt") 



# 2. Subgroup Analysis
PSM_High <- P4P_Clean_Binary %>% subset(PSM_C=="High")
PSM_Low <- P4P_Clean_Binary %>% subset(PSM_C=="Low")
Eff_High <- P4P_Clean_Binary %>% subset(Efficacy_C=="High")
Eff_Low <- P4P_Clean_Binary %>% subset(Efficacy_C=="Low")


# 2.1 Marginal Means
# 2.1.1 PSM
PSM_tp_High <- lm(Chosen_Job ~ 0 + `Total pay binary`, data=PSM_High)
PSM_tp_High.cl <- cluster.vcov(PSM_tp_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_tp_H <- broom::tidy(coeftest(PSM_tp_High, PSM_tp_High.cl))
PSM_tp_Low <- lm(Chosen_Job ~ 0 + `Total pay binary`, data=PSM_Low)
PSM_tp_Low.cl <- cluster.vcov(PSM_tp_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_tp_L <- broom::tidy(coeftest(PSM_tp_Low, PSM_tp_Low.cl))

PSM_pb_High <- lm(Chosen_Job ~ 0 + `Performance bonuses binary`, data=PSM_High)
PSM_pb_High.cl <- cluster.vcov(PSM_pb_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_pb_H <- broom::tidy(coeftest(PSM_pb_High, PSM_pb_High.cl))
PSM_pb_Low <- lm(Chosen_Job ~ 0 + `Performance bonuses binary`, data=PSM_Low)
PSM_pb_Low.cl <- cluster.vcov(PSM_pb_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_pb_L <- broom::tidy(coeftest(PSM_pb_Low, PSM_pb_Low.cl))

PSM_ge_High <- lm(Chosen_Job ~ 0 + `Goal based evaluation`, data=PSM_High)
PSM_ge_High.cl <- cluster.vcov(PSM_ge_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_ge_H <- broom::tidy(coeftest(PSM_ge_High, PSM_ge_High.cl))
PSM_ge_Low <- lm(Chosen_Job ~ 0 + `Goal based evaluation`, data=PSM_Low)
PSM_ge_Low.cl <- cluster.vcov(PSM_ge_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_ge_L <- broom::tidy(coeftest(PSM_ge_Low, PSM_ge_Low.cl))

PSM_cci_High <- lm(Chosen_Job ~ 0 + `Current community involvement binary`, data=PSM_High)
PSM_cci_High.cl <- cluster.vcov(PSM_cci_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_cci_H <- broom::tidy(coeftest(PSM_cci_High, PSM_cci_High.cl))
PSM_cci_Low <- lm(Chosen_Job ~ 0 + `Current community involvement binary`, data=PSM_Low)
PSM_cci_Low.cl <- cluster.vcov(PSM_cci_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_cci_L <- broom::tidy(coeftest(PSM_cci_Low, PSM_cci_Low.cl))

PSM_ci_High <- lm(Chosen_Job ~ 0 + `Community income binary`, data=PSM_High)
PSM_ci_High.cl <- cluster.vcov(PSM_ci_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_ci_H <- broom::tidy(coeftest(PSM_ci_High, PSM_ci_High.cl))
PSM_ci_Low <- lm(Chosen_Job ~ 0 + `Community income binary`, data=PSM_Low)
PSM_ci_Low.cl <- cluster.vcov(PSM_ci_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_ci_L <- broom::tidy(coeftest(PSM_ci_Low, PSM_ci_Low.cl))

PSM_cd_High <- lm(Chosen_Job ~ 0 + `Community demographics binary`, data=PSM_High)
PSM_cd_High.cl <- cluster.vcov(PSM_cd_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_cd_H <- broom::tidy(coeftest(PSM_cd_High, PSM_cd_High.cl))
PSM_cd_Low <- lm(Chosen_Job ~ 0 + `Community demographics binary`, data=PSM_Low)
PSM_cd_Low.cl <- cluster.vcov(PSM_cd_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_cd_L <- broom::tidy(coeftest(PSM_cd_Low, PSM_cd_Low.cl))

PSM_ow_High <- lm(Chosen_Job ~ 0 + `Overtime work binary`, data=PSM_High)
PSM_ow_High.cl <- cluster.vcov(PSM_ow_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_ow_H <- broom::tidy(coeftest(PSM_ow_High, PSM_ow_High.cl))
PSM_ow_Low <- lm(Chosen_Job ~ 0 + `Overtime work binary`, data=PSM_Low)
PSM_ow_Low.cl <- cluster.vcov(PSM_ow_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_ow_L <- broom::tidy(coeftest(PSM_ow_Low, PSM_ow_Low.cl))

PSM_kt_High <- lm(Chosen_Job ~ 0 + `Key job task binary`, data=PSM_High)
PSM_kt_High.cl <- cluster.vcov(PSM_kt_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_kt_H <- broom::tidy(coeftest(PSM_kt_High, PSM_kt_High.cl))
PSM_kt_Low <- lm(Chosen_Job ~ 0 + `Key job task binary`, data=PSM_Low)
PSM_kt_Low.cl <- cluster.vcov(PSM_kt_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_kt_L <- broom::tidy(coeftest(PSM_kt_Low, PSM_kt_Low.cl))

PSM_MM_H <- rbind(PSM_tp_H,PSM_pb_H,PSM_ge_H,PSM_cci_H,
                  PSM_ci_H,PSM_cd_H,PSM_ow_H,PSM_kt_H) %>%
  mutate(Group="High")

PSM_MM_L <- rbind(PSM_tp_L,PSM_pb_L,PSM_ge_L,PSM_cci_L,
                  PSM_ci_L,PSM_cd_L,PSM_ow_L,PSM_kt_L) %>%
  mutate(Group="Low")

PSM_MM <- rbind(PSM_MM_H,PSM_MM_L) %>%
  mutate(upper=estimate+qnorm(0.975)*std.error, lower=estimate-qnorm(0.975)*std.error) %>%
  mutate(Group=as.factor(Group))


# 2.1.2 Efficacy
Eff_tp_High <- lm(Chosen_Job ~ 0 + `Total pay binary`, data=Eff_High)
Eff_tp_High.cl <- cluster.vcov(Eff_tp_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_tp_H <- broom::tidy(coeftest(Eff_tp_High, Eff_tp_High.cl))
Eff_tp_Low <- lm(Chosen_Job ~ 0 + `Total pay binary`, data=Eff_Low)
Eff_tp_Low.cl <- cluster.vcov(Eff_tp_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_tp_L <- broom::tidy(coeftest(Eff_tp_Low, Eff_tp_Low.cl))

Eff_pb_High <- lm(Chosen_Job ~ 0 + `Performance bonuses binary`, data=Eff_High)
Eff_pb_High.cl <- cluster.vcov(Eff_pb_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_pb_H <- broom::tidy(coeftest(Eff_pb_High, Eff_pb_High.cl))
Eff_pb_Low <- lm(Chosen_Job ~ 0 + `Performance bonuses binary`, data=Eff_Low)
Eff_pb_Low.cl <- cluster.vcov(Eff_pb_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_pb_L <- broom::tidy(coeftest(Eff_pb_Low, Eff_pb_Low.cl))

Eff_ge_High <- lm(Chosen_Job ~ 0 + `Goal based evaluation`, data=Eff_High)
Eff_ge_High.cl <- cluster.vcov(Eff_ge_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_ge_H <- broom::tidy(coeftest(Eff_ge_High, Eff_ge_High.cl))
Eff_ge_Low <- lm(Chosen_Job ~ 0 + `Goal based evaluation`, data=Eff_Low)
Eff_ge_Low.cl <- cluster.vcov(Eff_ge_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_ge_L <- broom::tidy(coeftest(Eff_ge_Low, Eff_ge_Low.cl))

Eff_cci_High <- lm(Chosen_Job ~ 0 + `Current community involvement binary`, data=Eff_High)
Eff_cci_High.cl <- cluster.vcov(Eff_cci_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_cci_H <- broom::tidy(coeftest(Eff_cci_High, Eff_cci_High.cl))
Eff_cci_Low <- lm(Chosen_Job ~ 0 + `Current community involvement binary`, data=Eff_Low)
Eff_cci_Low.cl <- cluster.vcov(Eff_cci_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_cci_L <- broom::tidy(coeftest(Eff_cci_Low, Eff_cci_Low.cl))

Eff_ci_High <- lm(Chosen_Job ~ 0 + `Community income binary`, data=Eff_High)
Eff_ci_High.cl <- cluster.vcov(Eff_ci_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_ci_H <- broom::tidy(coeftest(Eff_ci_High, Eff_ci_High.cl))
Eff_ci_Low <- lm(Chosen_Job ~ 0 + `Community income binary`, data=Eff_Low)
Eff_ci_Low.cl <- cluster.vcov(Eff_ci_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_ci_L <- broom::tidy(coeftest(Eff_ci_Low, Eff_ci_Low.cl))

Eff_cd_High <- lm(Chosen_Job ~ 0 + `Community demographics binary`, data=Eff_High)
Eff_cd_High.cl <- cluster.vcov(Eff_cd_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_cd_H <- broom::tidy(coeftest(Eff_cd_High, Eff_cd_High.cl))
Eff_cd_Low <- lm(Chosen_Job ~ 0 + `Community demographics binary`, data=Eff_Low)
Eff_cd_Low.cl <- cluster.vcov(Eff_cd_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_cd_L <- broom::tidy(coeftest(Eff_cd_Low, Eff_cd_Low.cl))

Eff_ow_High <- lm(Chosen_Job ~ 0 + `Overtime work binary`, data=Eff_High)
Eff_ow_High.cl <- cluster.vcov(Eff_ow_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_ow_H <- broom::tidy(coeftest(Eff_ow_High, Eff_ow_High.cl))
Eff_ow_Low <- lm(Chosen_Job ~ 0 + `Overtime work binary`, data=Eff_Low)
Eff_ow_Low.cl <- cluster.vcov(Eff_ow_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_ow_L <- broom::tidy(coeftest(Eff_ow_Low, Eff_ow_Low.cl))

Eff_kt_High <- lm(Chosen_Job ~ 0 + `Key job task binary`, data=Eff_High)
Eff_kt_High.cl <- cluster.vcov(Eff_kt_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_kt_H <- broom::tidy(coeftest(Eff_kt_High, Eff_kt_High.cl))
Eff_kt_Low <- lm(Chosen_Job ~ 0 + `Key job task binary`, data=Eff_Low)
Eff_kt_Low.cl <- cluster.vcov(Eff_kt_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_kt_L <- broom::tidy(coeftest(Eff_kt_Low, Eff_kt_Low.cl))

Eff_MM_H <- rbind(Eff_tp_H,Eff_pb_H,Eff_ge_H,Eff_cci_H,
                  Eff_ci_H,Eff_cd_H,Eff_ow_H,Eff_kt_H) %>%
  mutate(Group="High")

Eff_MM_L <- rbind(Eff_tp_L,Eff_pb_L,Eff_ge_L,Eff_cci_L,
                  Eff_ci_L,Eff_cd_L,Eff_ow_L,Eff_kt_L) %>%
  mutate(Group="Low")

Eff_MM <- rbind(Eff_MM_H,Eff_MM_L) %>%
  mutate(upper=estimate+qnorm(0.975)*std.error, lower=estimate-qnorm(0.975)*std.error) %>%
  mutate(Group=as.factor(Group))


# 2.1.3 Plot (Figure 3)
PSM_MM.Match <- data.frame(trimws(gsub(".*[`]","",PSM_MM$term),"both"))
names(PSM_MM.Match) <- "Value"
PSM_MM.Match1 <- cbind(PSM_MM,PSM_MM.Match)

PSM.MM.table <- merge(x=Base.Cat, y=PSM_MM.Match1,
                      by='Value', all.x=TRUE)
PSM.MM.table <- within(PSM.MM.table, Group <- relevel(Group, ref = "Low"))

PSM.MM.table$Value <- factor(PSM.MM.table$Value,
                             levels=c("Total pay:", "Slightly below average","About average and Slightly above average",
                                      "Performance bonuses binary:", "No (fixed salary)","Yes (bonuses)",
                                      "Goal based evaluation:","No (supervisor based)","Yes (goal based)",
                                      "Current community involvement:","Rare participation","Moderate and Frequent participation",
                                      "Community income:","Low income","Average and High income",
                                      "Community demographics:","Mostly white","Non-White",
                                      "Overtime work:","Never required","Occasionally and Frequently required",
                                      "Key job task:","Analysis identifying community needs","Interaction with community and organizational members"))


Eff_MM.Match <- data.frame(trimws(gsub(".*[`]","",Eff_MM$term),"both"))
names(Eff_MM.Match) <- "Value"
Eff_MM.Match1 <- cbind(Eff_MM,Eff_MM.Match)

Eff.MM.table <- merge(x=Base.Cat, y=Eff_MM.Match1,
                      by='Value', all.x=TRUE)
Eff.MM.table <- within(Eff.MM.table, Group <- relevel(Group, ref = "Low"))

Eff.MM.table$Value <- factor(Eff.MM.table$Value,
                             levels=c("Total pay:", "Slightly below average","About average and Slightly above average",
                                      "Performance bonuses binary:", "No (fixed salary)","Yes (bonuses)",
                                      "Goal based evaluation:","No (supervisor based)","Yes (goal based)",
                                      "Current community involvement:","Rare participation","Moderate and Frequent participation",
                                      "Community income:","Low income","Average and High income",
                                      "Community demographics:","Mostly white","Non-White",
                                      "Overtime work:","Never required","Occasionally and Frequently required",
                                      "Key job task:","Analysis identifying community needs","Interaction with community and organizational members"))


PSM.MM.table.match <- PSM.MM.table %>%
  subset(Value!="No (fixed salary)" & Value!="Yes (bonuses)" & Value!="No (supervisor based)" & Value!="Yes (goal based)")
Eff.MM.table.match <- Eff.MM.table %>%
  subset(Value=="No (fixed salary)" | Value=="Yes (bonuses)" | Value=="No (supervisor based)" | Value=="Yes (goal based)")
PSM.Eff.MM.table <- rbind(PSM.MM.table.match,Eff.MM.table.match) %>%
  mutate(Value.fix=ifelse(Value=="Total pay:","Total pay: PSM",
                   ifelse(Value=="Performance bonuses binary:", "Performance bonuses binary: Self-Efficacy",
                   ifelse(Value=="Goal based evaluation:", "Goal based evaluation: Self-Efficacy",
                   ifelse(Value=="Current community involvement:","Current community involvement: PSM",
                   ifelse(Value=="Community income:","Community income: PSM",
                   ifelse(Value=="Community demographics:","Community demographics: PSM",
                   ifelse(Value=="Overtime work:","Overtime work: PSM",
                   ifelse(Value=="Key job task:","Key job task: PSM", as.character(Value))))))))))
PSM.Eff.MM.table$Value.fix <- factor(PSM.Eff.MM.table$Value.fix,
                                     levels=c("Total pay: PSM", "Slightly below average","About average and Slightly above average",
                                              "Performance bonuses binary: Self-Efficacy", "No (fixed salary)","Yes (bonuses)",
                                              "Goal based evaluation: Self-Efficacy","No (supervisor based)","Yes (goal based)",
                                              "Current community involvement: PSM","Rare participation","Moderate and Frequent participation",
                                              "Community income: PSM","Low income","Average and High income",
                                              "Community demographics: PSM","Mostly white","Non-White",
                                              "Overtime work: PSM","Never required","Occasionally and Frequently required",
                                              "Key job task: PSM","Analysis identifying community needs","Interaction with community and organizational members"))

ggplot(data=PSM.Eff.MM.table,
       aes(estimate, Value.fix, color=Group))+
  geom_point(position = position_dodge(width = 0.6), size = 1.5) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), 
                 position = position_dodge(width = 0.6), size = .3) +
  geom_vline(aes(xintercept = 0.5), size = .25, linetype = "dashed") +
  scale_y_discrete(limits=rev) +
  theme(axis.text.y = element_text(hjust=0)) +
  theme(axis.text.y = element_text(face = 
                                     c("plain","plain","bold","plain","plain","bold","plain","plain",
                                       "bold","plain","plain","bold","plain","plain","bold","plain","plain","bold","plain","plain","bold",
                                       "plain","plain","bold"))) +
  ylab("") +
  xlab("Marginal Mean") +
  theme(axis.text.y = element_text(size=12, colour = "black")) +
  theme(text = element_text(size=12, colour = "black")) + 
  scale_colour_discrete(na.translate = F,limits=c("High","Low")) +
  theme(legend.position = "bottom") + 
  theme(legend.title=element_blank())


# 2.1.4 Differences in marginal means (Appendix Figure E2)
PSM_MM_Diff <- merge(x=PSM_MM_H, y=PSM_MM_L, by=c('term'), all.x=TRUE)
PSM_MM_Diff1 <- PSM_MM_Diff %>%
  mutate(estimate_diff=estimate.x-estimate.y) %>%
  mutate(std.error_diff=sqrt((std.error.x^2)+(std.error.y^2))) %>%
  mutate(upper=estimate_diff+qnorm(0.975)*std.error_diff, lower=estimate_diff-qnorm(0.975)*std.error_diff) 
PSM_MM_Diff.Match <- data.frame(trimws(gsub(".*[`]","",PSM_MM_Diff1$term),"both"))
names(PSM_MM_Diff.Match) <- "Value"
PSM_MM_Diff.Match1 <- cbind(PSM_MM_Diff1,PSM_MM_Diff.Match)
PSM_MM_Diff.table <- merge(x=Base.Cat, y=PSM_MM_Diff.Match1,
                           by='Value', all.x=TRUE)

Eff_MM_Diff <- merge(x=Eff_MM_H, y=Eff_MM_L, by=c('term'), all.x=TRUE)
Eff_MM_Diff1 <- Eff_MM_Diff %>%
  mutate(estimate_diff=estimate.x-estimate.y) %>%
  mutate(std.error_diff=sqrt((std.error.x^2)+(std.error.y^2))) %>%
  mutate(upper=estimate_diff+qnorm(0.975)*std.error_diff, lower=estimate_diff-qnorm(0.975)*std.error_diff) 
Eff_MM_Diff.Match <- data.frame(trimws(gsub(".*[`]","",Eff_MM_Diff1$term),"both"))
names(Eff_MM_Diff.Match) <- "Value"
Eff_MM_Diff.Match1 <- cbind(Eff_MM_Diff1,Eff_MM_Diff.Match)
Eff_MM_Diff.table <- merge(x=Base.Cat, y=Eff_MM_Diff.Match1,
                           by='Value', all.x=TRUE)

PSM_MM_Diff.match <- PSM_MM_Diff.table %>%
  subset(Value!="No (fixed salary)" & Value!="Yes (bonuses)" & Value!="No (supervisor based)" & Value!="Yes (goal based)")
Eff_MM_Diff.match <- Eff_MM_Diff.table %>%
  subset(Value=="No (fixed salary)" | Value=="Yes (bonuses)" | Value=="No (supervisor based)" | Value=="Yes (goal based)")
PSM.Eff.MM.Diff.table <- rbind(PSM_MM_Diff.match,Eff_MM_Diff.match) %>%
  mutate(Value.fix=ifelse(Value=="Total pay:","Total pay: PSM",
                   ifelse(Value=="Performance bonuses binary:", "Performance bonuses binary: Self-Efficacy",
                   ifelse(Value=="Goal based evaluation:", "Goal based evaluation: Self-Efficacy",
                   ifelse(Value=="Current community involvement:","Current community involvement: PSM",
                   ifelse(Value=="Community income:","Community income: PSM",
                   ifelse(Value=="Community demographics:","Community demographics: PSM",
                   ifelse(Value=="Overtime work:","Overtime work: PSM",
                   ifelse(Value=="Key job task:","Key job task: PSM", as.character(Value))))))))))
PSM.Eff.MM.Diff.table$Value.fix <- factor(PSM.Eff.MM.Diff.table$Value.fix,
                                     levels=c("Total pay: PSM", "Slightly below average","About average and Slightly above average",
                                              "Performance bonuses binary: Self-Efficacy", "No (fixed salary)","Yes (bonuses)",
                                              "Goal based evaluation: Self-Efficacy","No (supervisor based)","Yes (goal based)",
                                              "Current community involvement: PSM","Rare participation","Moderate and Frequent participation",
                                              "Community income: PSM","Low income","Average and High income",
                                              "Community demographics: PSM","Mostly white","Non-White",
                                              "Overtime work: PSM","Never required","Occasionally and Frequently required",
                                              "Key job task: PSM","Analysis identifying community needs","Interaction with community and organizational members"))

ggplot(PSM.Eff.MM.Diff.table, aes(x = estimate_diff, y = Value.fix)) + 
  geom_vline(aes(xintercept = 0.0), size = .25, linetype = "dashed") + 
  geom_errorbarh(data = PSM.Eff.MM.Diff.table, aes(xmax = upper, xmin = lower), size = .3, height = .6, position = position_nudge(y = 0), color="black") +
  geom_point(data = PSM.Eff.MM.Diff.table, shape=19, size = 1.5, color='black') +
  scale_y_discrete(limits=rev) +
  theme(axis.text.y = element_text(hjust=0)) +
  theme(axis.text.y = element_text(face = 
                                     c("plain","plain","bold","plain","plain","bold","plain","plain",
                                       "bold","plain","plain","bold","plain","plain","bold","plain","plain","bold","plain","plain","bold",
                                       "plain","plain","bold"))) +
  ylab("") +
  xlab("Difference in Marginal Means") +
  theme(axis.text.y = element_text(size=12, colour = "black")) +
  theme(text = element_text(size=12, colour = "black"))



# 2.1.5 Subgroup MM (all levels)
PSM_tp_High <- lm(Chosen_Job ~ 0 + `Total pay`, data=PSM_High)
PSM_tp_High.cl <- cluster.vcov(PSM_tp_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_tp_H <- broom::tidy(coeftest(PSM_tp_High, PSM_tp_High.cl))
PSM_tp_Low <- lm(Chosen_Job ~ 0 + `Total pay`, data=PSM_Low)
PSM_tp_Low.cl <- cluster.vcov(PSM_tp_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_tp_L <- broom::tidy(coeftest(PSM_tp_Low, PSM_tp_Low.cl))

PSM_pb_High <- lm(Chosen_Job ~ 0 + `Performance bonuses`, data=PSM_High)
PSM_pb_High.cl <- cluster.vcov(PSM_pb_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_pb_H <- broom::tidy(coeftest(PSM_pb_High, PSM_pb_High.cl))
PSM_pb_Low <- lm(Chosen_Job ~ 0 + `Performance bonuses`, data=PSM_Low)
PSM_pb_Low.cl <- cluster.vcov(PSM_pb_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_pb_L <- broom::tidy(coeftest(PSM_pb_Low, PSM_pb_Low.cl))

PSM_ge_High <- lm(Chosen_Job ~ 0 + `Job performance evaluation`, data=PSM_High)
PSM_ge_High.cl <- cluster.vcov(PSM_ge_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_ge_H <- broom::tidy(coeftest(PSM_ge_High, PSM_ge_High.cl))
PSM_ge_Low <- lm(Chosen_Job ~ 0 + `Job performance evaluation`, data=PSM_Low)
PSM_ge_Low.cl <- cluster.vcov(PSM_ge_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_ge_L <- broom::tidy(coeftest(PSM_ge_Low, PSM_ge_Low.cl))

PSM_cci_High <- lm(Chosen_Job ~ 0 + `Current community involvement`, data=PSM_High)
PSM_cci_High.cl <- cluster.vcov(PSM_cci_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_cci_H <- broom::tidy(coeftest(PSM_cci_High, PSM_cci_High.cl))
PSM_cci_Low <- lm(Chosen_Job ~ 0 + `Current community involvement`, data=PSM_Low)
PSM_cci_Low.cl <- cluster.vcov(PSM_cci_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_cci_L <- broom::tidy(coeftest(PSM_cci_Low, PSM_cci_Low.cl))

PSM_ci_High <- lm(Chosen_Job ~ 0 + `Community income`, data=PSM_High)
PSM_ci_High.cl <- cluster.vcov(PSM_ci_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_ci_H <- broom::tidy(coeftest(PSM_ci_High, PSM_ci_High.cl))
PSM_ci_Low <- lm(Chosen_Job ~ 0 + `Community income`, data=PSM_Low)
PSM_ci_Low.cl <- cluster.vcov(PSM_ci_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_ci_L <- broom::tidy(coeftest(PSM_ci_Low, PSM_ci_Low.cl))

PSM_cd_High <- lm(Chosen_Job ~ 0 + `Community demographics`, data=PSM_High)
PSM_cd_High.cl <- cluster.vcov(PSM_cd_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_cd_H <- broom::tidy(coeftest(PSM_cd_High, PSM_cd_High.cl))
PSM_cd_Low <- lm(Chosen_Job ~ 0 + `Community demographics`, data=PSM_Low)
PSM_cd_Low.cl <- cluster.vcov(PSM_cd_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_cd_L <- broom::tidy(coeftest(PSM_cd_Low, PSM_cd_Low.cl))

PSM_ow_High <- lm(Chosen_Job ~ 0 + `Overtime work`, data=PSM_High)
PSM_ow_High.cl <- cluster.vcov(PSM_ow_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_ow_H <- broom::tidy(coeftest(PSM_ow_High, PSM_ow_High.cl))
PSM_ow_Low <- lm(Chosen_Job ~ 0 + `Overtime work`, data=PSM_Low)
PSM_ow_Low.cl <- cluster.vcov(PSM_ow_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_ow_L <- broom::tidy(coeftest(PSM_ow_Low, PSM_ow_Low.cl))

PSM_kt_High <- lm(Chosen_Job ~ 0 + `Key job task`, data=PSM_High)
PSM_kt_High.cl <- cluster.vcov(PSM_kt_High, cbind(PSM_High$ID,PSM_High$Treatment.ID))
PSM_kt_H <- broom::tidy(coeftest(PSM_kt_High, PSM_kt_High.cl))
PSM_kt_Low <- lm(Chosen_Job ~ 0 + `Key job task`, data=PSM_Low)
PSM_kt_Low.cl <- cluster.vcov(PSM_kt_Low, cbind(PSM_Low$ID,PSM_Low$Treatment.ID))
PSM_kt_L <- broom::tidy(coeftest(PSM_kt_Low, PSM_kt_Low.cl))

PSM_MM_H <- rbind(PSM_tp_H,PSM_pb_H,PSM_ge_H,PSM_cci_H,
                PSM_ci_H,PSM_cd_H,PSM_ow_H,PSM_kt_H) %>%
  mutate(Group="High")

PSM_MM_L <- rbind(PSM_tp_L,PSM_pb_L,PSM_ge_L,PSM_cci_L,
                PSM_ci_L,PSM_cd_L,PSM_ow_L,PSM_kt_L) %>%
  mutate(Group="Low")

PSM_MM <- rbind(PSM_MM_H,PSM_MM_L) %>%
  mutate(upper=estimate+qnorm(0.975)*std.error, lower=estimate-qnorm(0.975)*std.error) %>%
  mutate(Group=as.factor(Group))


Eff_tp_High <- lm(Chosen_Job ~ 0 + `Total pay`, data=Eff_High)
Eff_tp_High.cl <- cluster.vcov(Eff_tp_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_tp_H <- broom::tidy(coeftest(Eff_tp_High, Eff_tp_High.cl))
Eff_tp_Low <- lm(Chosen_Job ~ 0 + `Total pay`, data=Eff_Low)
Eff_tp_Low.cl <- cluster.vcov(Eff_tp_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_tp_L <- broom::tidy(coeftest(Eff_tp_Low, Eff_tp_Low.cl))

Eff_pb_High <- lm(Chosen_Job ~ 0 + `Performance bonuses`, data=Eff_High)
Eff_pb_High.cl <- cluster.vcov(Eff_pb_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_pb_H <- broom::tidy(coeftest(Eff_pb_High, Eff_pb_High.cl))
Eff_pb_Low <- lm(Chosen_Job ~ 0 + `Performance bonuses`, data=Eff_Low)
Eff_pb_Low.cl <- cluster.vcov(Eff_pb_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_pb_L <- broom::tidy(coeftest(Eff_pb_Low, Eff_pb_Low.cl))

Eff_ge_High <- lm(Chosen_Job ~ 0 + `Job performance evaluation`, data=Eff_High)
Eff_ge_High.cl <- cluster.vcov(Eff_ge_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_ge_H <- broom::tidy(coeftest(Eff_ge_High, Eff_ge_High.cl))
Eff_ge_Low <- lm(Chosen_Job ~ 0 + `Job performance evaluation`, data=Eff_Low)
Eff_ge_Low.cl <- cluster.vcov(Eff_ge_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_ge_L <- broom::tidy(coeftest(Eff_ge_Low, Eff_ge_Low.cl))

Eff_cci_High <- lm(Chosen_Job ~ 0 + `Current community involvement`, data=Eff_High)
Eff_cci_High.cl <- cluster.vcov(Eff_cci_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_cci_H <- broom::tidy(coeftest(Eff_cci_High, Eff_cci_High.cl))
Eff_cci_Low <- lm(Chosen_Job ~ 0 + `Current community involvement`, data=Eff_Low)
Eff_cci_Low.cl <- cluster.vcov(Eff_cci_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_cci_L <- broom::tidy(coeftest(Eff_cci_Low, Eff_cci_Low.cl))

Eff_ci_High <- lm(Chosen_Job ~ 0 + `Community income`, data=Eff_High)
Eff_ci_High.cl <- cluster.vcov(Eff_ci_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_ci_H <- broom::tidy(coeftest(Eff_ci_High, Eff_ci_High.cl))
Eff_ci_Low <- lm(Chosen_Job ~ 0 + `Community income`, data=Eff_Low)
Eff_ci_Low.cl <- cluster.vcov(Eff_ci_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_ci_L <- broom::tidy(coeftest(Eff_ci_Low, Eff_ci_Low.cl))

Eff_cd_High <- lm(Chosen_Job ~ 0 + `Community demographics`, data=Eff_High)
Eff_cd_High.cl <- cluster.vcov(Eff_cd_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_cd_H <- broom::tidy(coeftest(Eff_cd_High, Eff_cd_High.cl))
Eff_cd_Low <- lm(Chosen_Job ~ 0 + `Community demographics`, data=Eff_Low)
Eff_cd_Low.cl <- cluster.vcov(Eff_cd_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_cd_L <- broom::tidy(coeftest(Eff_cd_Low, Eff_cd_Low.cl))

Eff_ow_High <- lm(Chosen_Job ~ 0 + `Overtime work`, data=Eff_High)
Eff_ow_High.cl <- cluster.vcov(Eff_ow_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_ow_H <- broom::tidy(coeftest(Eff_ow_High, Eff_ow_High.cl))
Eff_ow_Low <- lm(Chosen_Job ~ 0 + `Overtime work`, data=Eff_Low)
Eff_ow_Low.cl <- cluster.vcov(Eff_ow_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_ow_L <- broom::tidy(coeftest(Eff_ow_Low, Eff_ow_Low.cl))

Eff_kt_High <- lm(Chosen_Job ~ 0 + `Key job task`, data=Eff_High)
Eff_kt_High.cl <- cluster.vcov(Eff_kt_High, cbind(Eff_High$ID,Eff_High$Treatment.ID))
Eff_kt_H <- broom::tidy(coeftest(Eff_kt_High, Eff_kt_High.cl))
Eff_kt_Low <- lm(Chosen_Job ~ 0 + `Key job task`, data=Eff_Low)
Eff_kt_Low.cl <- cluster.vcov(Eff_kt_Low, cbind(Eff_Low$ID,Eff_Low$Treatment.ID))
Eff_kt_L <- broom::tidy(coeftest(Eff_kt_Low, Eff_kt_Low.cl))

Eff_MM_H <- rbind(Eff_tp_H,Eff_pb_H,Eff_ge_H,Eff_cci_H,
                  Eff_ci_H,Eff_cd_H,Eff_ow_H,Eff_kt_H) %>%
  mutate(Group="High")

Eff_MM_L <- rbind(Eff_tp_L,Eff_pb_L,Eff_ge_L,Eff_cci_L,
                  Eff_ci_L,Eff_cd_L,Eff_ow_L,Eff_kt_L) %>%
  mutate(Group="Low")

Eff_MM <- rbind(Eff_MM_H,Eff_MM_L) %>%
  mutate(upper=estimate+qnorm(0.975)*std.error, lower=estimate-qnorm(0.975)*std.error) %>%
  mutate(Group=as.factor(Group))


# Plot (Appendix Figure F1)
Base.Cat <- data.frame(
  Value=c("Total pay:", "Slightly below average","About average","Slightly above average",
          "Performance bonuses:", "No performance bonuses; fixed salary","A small part of your potential pay (5%)","A moderate part of your potential pay (10%)","A large part of your potential pay (20%)",
          "Job performance evaluation:","A supervisor evaluation of your work","Attendance numbers for Project HOPE events","Satisfaction surveys of Project HOPE event participants","Changes in community crime, poverty, and blight",
          "Current community involvement:","Rare participation","Moderate participation","Frequent participation",
          "Community income:","Low income","Average income","High income",
          "Community demographics:","Mostly white","Mostly African American","Mostly Hispanic","Multiracial",
          "Overtime work:","Never required","Occasionally required","Frequently required",
          "Key job task:","Analysis identifying community needs","Coordination with community groups and organizations",
          "Direct interaction with community residents","Teamwork with peers and supervisors"))

PSM_MM.Match <- data.frame(trimws(gsub(".*[`]","",PSM_MM$term),"both"))
names(PSM_MM.Match) <- "Value"
PSM_MM.Match1 <- cbind(PSM_MM,PSM_MM.Match)

PSM.MM.table <- merge(x=Base.Cat, y=PSM_MM.Match1,
                       by='Value', all.x=TRUE)
PSM.MM.table <- within(PSM.MM.table, Group <- relevel(Group, ref = "Low"))


Eff_MM.Match <- data.frame(trimws(gsub(".*[`]","",Eff_MM$term),"both"))
names(Eff_MM.Match) <- "Value"
Eff_MM.Match1 <- cbind(Eff_MM,Eff_MM.Match)

Eff.MM.table <- merge(x=Base.Cat, y=Eff_MM.Match1,
                      by='Value', all.x=TRUE)
Eff.MM.table <- within(Eff.MM.table, Group <- relevel(Group, ref = "Low"))


PSM.MM.table.match <- PSM.MM.table %>%
  subset(Value!="No performance bonuses; fixed salary" & Value!="A small part of your potential pay (5%)" & Value!="A moderate part of your potential pay (10%)" & Value!="A large part of your potential pay (20%)" & Value!="A supervisor evaluation of your work" & Value!="Attendance numbers for Project HOPE events" & Value!="Satisfaction surveys of Project HOPE event participants" & Value!="Changes in community crime, poverty, and blight")
Eff.MM.table.match <- Eff.MM.table %>%
  subset(Value=="No performance bonuses; fixed salary" | Value=="A small part of your potential pay (5%)" | Value=="A moderate part of your potential pay (10%)" | Value=="A large part of your potential pay (20%)" | Value=="A supervisor evaluation of your work" | Value=="Attendance numbers for Project HOPE events" | Value=="Satisfaction surveys of Project HOPE event participants" | Value=="Changes in community crime, poverty, and blight")

PSM.Eff.MM.table <- rbind(PSM.MM.table.match,Eff.MM.table.match) %>%
  mutate(Value.fix=ifelse(Value=="Total pay:","Total pay: PSM",
                   ifelse(Value=="Performance bonuses:", "Performance bonuses: Self-Efficacy",
                   ifelse(Value=="Job performance evaluation:", "Job performance evaluation: Self-Efficacy",
                   ifelse(Value=="Current community involvement:","Current community involvement: PSM",
                   ifelse(Value=="Community income:","Community income: PSM",
                   ifelse(Value=="Community demographics:","Community demographics: PSM",
                   ifelse(Value=="Overtime work:","Overtime work: PSM",
                   ifelse(Value=="Key job task:","Key job task: PSM", as.character(Value))))))))))
PSM.Eff.MM.table$Value.fix <- factor(PSM.Eff.MM.table$Value.fix,
                             levels=c("Total pay: PSM", "Slightly below average","About average","Slightly above average",
                                      "Performance bonuses: Self-Efficacy","No performance bonuses; fixed salary", "A small part of your potential pay (5%)","A moderate part of your potential pay (10%)","A large part of your potential pay (20%)",
                                      "Job performance evaluation: Self-Efficacy","A supervisor evaluation of your work","Attendance numbers for Project HOPE events","Satisfaction surveys of Project HOPE event participants","Changes in community crime, poverty, and blight",
                                      "Current community involvement: PSM","Rare participation","Moderate participation","Frequent participation",
                                      "Community income: PSM","Low income","Average income","High income",
                                      "Community demographics: PSM","Mostly white","Mostly African American","Mostly Hispanic","Multiracial",
                                      "Overtime work: PSM","Never required","Occasionally required","Frequently required",
                                      "Key job task: PSM","Analysis identifying community needs","Coordination with community groups and organizations",
                                      "Direct interaction with community residents","Teamwork with peers and supervisors"))

ggplot(data=PSM.Eff.MM.table,
       aes(estimate, Value.fix, color=Group))+
  geom_point(position = position_dodge(width = 0.6), size = 1.5) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), 
                 position = position_dodge(width = 0.6), size = .3) +
  geom_vline(aes(xintercept = 0.5), size = .25, linetype = "dashed") +
  scale_y_discrete(limits=rev) +
  theme(axis.text.y = element_text(hjust=0)) +
  theme(axis.text.y = element_text(face = 
                                     c("plain","plain","plain","plain","bold","plain","plain","plain","bold","plain","plain","plain","plain","bold",
                                       "plain","plain","plain","bold","plain","plain","plain","bold","plain","plain","plain","plain","bold","plain","plain","plain","plain","bold",
                                       "plain","plain","plain","bold"))) +
  ylab("") +
  xlab("Marginal Mean") +
  theme(axis.text.y = element_text(size=12, colour = "black")) +
  theme(text = element_text(size=12, colour = "black")) + 
  scale_colour_discrete(na.translate = F,limits=c("High","Low")) +
  theme(legend.position = "bottom") + 
  theme(legend.title=element_blank())


# 2.2 AMCE 
# 2.2.1. Full Model (Appendix Table E1)
Full.AMCE <- lm(Chosen_Job ~ `Total pay binary`*PSM_C + `Performance bonuses binary`*Efficacy_C + `Goal based evaluation`*Efficacy_C +
                  `Current community involvement binary`*PSM_C + `Community income binary`*PSM_C + `Community demographics binary`*PSM_C +
                  `Overtime work binary`*PSM_C + `Key job task binary`*PSM_C, 
                data=P4P_Clean_Binary)
summary(Full.AMCE)
Full.AMCE.cl <- cluster.vcov(Full.AMCE, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
Full.stat <- coeftest(Full.AMCE, Full.AMCE.cl)
Twoway.se <- sqrt(diag(Full.AMCE.cl))

# BH Procedure
pval <- Full.stat[,4]
p.bh <- p.adjust(pval, "fdr")
p.bh05 <- ifelse(p.bh<=0.05, 1, 0)
p.bh10 <- ifelse(p.bh<=0.10, 1, 0)


Full.AMCE.Collapse <- lm(Chosen_Job ~ `Total pay`*PSM_C + `Performance bonuses`*Efficacy_C + `Job performance evaluation`*Efficacy_C +
                  `Current community involvement`*PSM_C + `Community income`*PSM_C + `Community demographics`*PSM_C +
                  `Overtime work`*PSM_C + `Key job task`*PSM_C, 
                data=P4P_Clean_Binary)
summary(Full.AMCE.Collapse)
Full.AMCE.Collapse.cl <- cluster.vcov(Full.AMCE.Collapse, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
Full.Collapse.stat <- coeftest(Full.AMCE.Collapse, Full.AMCE.Collapse.cl)
Twoway.se.Collapse <- sqrt(diag(Full.AMCE.Collapse.cl))

# BH Procedure
pval <- Full.Collapse.stat[,4]
p.bh <- p.adjust(pval, "fdr")
p.bh05 <- ifelse(p.bh<=0.05, 1, 0)
p.bh10 <- ifelse(p.bh<=0.10, 1, 0)


# Reporting
stargazer(Full.AMCE,Full.AMCE.Collapse,
          se=list(Twoway.se,Twoway.se.Collapse), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          order = c(1,2,3,5,6,7,8,10,24,25,26,11,27,28,12,29,30,13,31,32,33,14,34,35,15,36,37,38,
                    4,9,16,39,40,17,41,42,43,18,44,45,46,19,47,48,20,49,50,21,51,52,53,22,54,55,23,56,57,58),
          covariate.labels=c("Total pay binary: About average and Slightly above average",
                             "Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes",
                             "Performance bonuses: large (20%)",
                             "Performance bonuses: moderate (10%)",
                             "Performance bonuses: small (5%)",
                             "Goal based evaluation: Yes",
                             "Goal based evaluation: Attendance",
                             "Goal based evaluation: Changes",
                             "Goal based evaluation: Satisfaction",
                             "Current community involvement binary: Moderate and Frequent participation",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income binary: Average and High income",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics binary: Non-White",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work binary: Occasionally and Frequently required",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task binary: Interaction with community and organizational members",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "PSM (A)",
                             "Self-Efficacy (B)",
                             "Total pay binary: About average and Slightly above average X (A)",
                             "Total pay: About average X (A)","Total pay: Slightly above average X (A)",
                             "Performance bonuses binary: Yes X (B)",
                             "Performance bonuses: large (20%)  X (B)",
                             "Performance bonuses: moderate (10%) X (B)",
                             "Performance bonuses: small (5%) X (B)",
                             "Goal based evaluation: Yes X (B)",
                             "Goal based evaluation: Attendance X (B)",
                             "Goal based evaluation: Changes X (B)",
                             "Goal based evaluation: Satisfaction X (B)",
                             "Current community involvement binary: Moderate and Frequent participation X (A)",
                             "Current community involvement: Frequent participation X (A)",
                             "Current community involvement: Moderate participation X (A)",
                             "Community income binary: Average and High income X (A)",
                             "Community income: Average income X (A)", "Community income: High income X (A)",
                             "Community demographics binary: Non-White X (A)",
                             "Community demographics: Mostly African American X (A)",
                             "Community demographics: Mostly Hispanic X (A)",
                             "Community demographics: Multiracial X (A)",
                             "Overtime work binary: Occasionally and Frequently required X (A)",
                             "Overtime work: Frequently required X (A)","Overtime work: Occasionally required X (A)",
                             "Key job task binary: Interaction with community and organizational members X (A)",
                             "Key job task: Coordination with community groups X (A)",
                             "Key job task: Direct interaction with community residents X (A)",
                             "Key job task: Teamwork with peers and supervisors X (A)"),
          out="Appendix_Table_E1.txt") 


# Plot (Appendix Figure E1)
Interaction.Full.PSM <- margins(Full.AMCE, 
                                vcov = cluster.vcov(Full.AMCE, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID)), 
                                variables = c("Total pay binary","Current community involvement binary","Community income binary",
                                              "Community demographics binary","Overtime work binary","Key job task binary"),
                                at=list(PSM_C=c("High","Low"))) %>%
  summary() %>% rename(Value=PSM_C) %>% 
  mutate(group=ifelse(factor=="Total pay binaryAbout average and Slightly above average" & Value=="High","About average and Slightly above average: High PSM",
               ifelse(factor=="Current community involvement binaryModerate and Frequent participation" & Value=="High", "Moderate and Frequent participation: High PSM",
               ifelse(factor=="Community income binaryAverage and High income" & Value=="High", "Average and High income: High PSM",
               ifelse(factor=="Community demographics binaryNon-White" & Value=="High", "Non-White: High PSM",
               ifelse(factor=="Key job task binaryInteraction with community and organizational members" & Value=="High","Interaction with community and organizational members: High PSM",
               ifelse(factor=="Overtime work binaryOccasionally and Frequently required" & Value=="High", "Occasionally and Frequently required: High PSM",
               ifelse(factor=="Total pay binaryAbout average and Slightly above average" & Value=="Low","About average and Slightly above average: Low PSM",
               ifelse(factor=="Current community involvement binaryModerate and Frequent participation" & Value=="Low", "Moderate and Frequent participation: Low PSM",
               ifelse(factor=="Community income binaryAverage and High income" & Value=="Low", "Average and High income: Low PSM",
               ifelse(factor=="Community demographics binaryNon-White" & Value=="Low", "Non-White: Low PSM",
               ifelse(factor=="Key job task binaryInteraction with community and organizational members" & Value=="Low","Interaction with community and organizational members: Low PSM",
               ifelse(factor=="Overtime work binaryOccasionally and Frequently required" & Value=="Low", "Occasionally and Frequently required: Low PSM","XXX"))))))))))))) %>%
  mutate(Int="PSM") %>%
  mutate(upper_90=AME+qnorm(0.95)*SE, lower_90=AME-qnorm(0.95)*SE)

Interaction.Full.PSE <- margins(Full.AMCE, 
                                vcov = cluster.vcov(Full.AMCE, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID)), 
                                variables = c("Performance bonuses binary","Goal based evaluation"),
                                at=list(Efficacy_C=c("High","Low"))) %>%
  summary() %>% rename(Value=Efficacy_C) %>% 
  mutate(group=ifelse(factor=="Performance bonuses binaryYes (bonuses)" & Value=="High","Yes (bonuses): High Self-Efficacy",
               ifelse(factor=="Goal based evaluationYes (goal based)" & Value=="High","Yes (goal based): High Self-Efficacy",
               ifelse(factor=="Performance bonuses binaryYes (bonuses)" & Value=="Low","Yes (bonuses): Low Self-Efficacy",
               ifelse(factor=="Goal based evaluationYes (goal based)" & Value=="Low","Yes (goal based): Low Self-Efficacy","XXX"))))) %>%
  mutate(Int="Self-efficacy") %>%
  mutate(upper_90=AME+qnorm(0.95)*SE, lower_90=AME-qnorm(0.95)*SE)

Interaction.Full <- rbind(Interaction.Full.PSM,Interaction.Full.PSE) 

Base.Factor <- data.frame(
  group=c("Total pay:", "(Baseline=Slightly below aveage)","About average and Slightly above average: High PSM","About average and Slightly above average: Low PSM",
          "Performance bonuses binary:", "(Baseline=No (Fixed salary))","Yes (bonuses): High Self-Efficacy","Yes (bonuses): Low Self-Efficacy",
          "Goal based evaluation:","(Baseline=No (Supervisor based)","Yes (goal based): High Self-Efficacy","Yes (goal based): Low Self-Efficacy",
          "Current community involvement:","(Baseline=Rare participation)","Moderate and Frequent participation: High PSM","Moderate and Frequent participation: Low PSM",
          "Community income:","(Baseline=Low income)","Average and High income: High PSM","Average and High income: Low PSM",
          "Community demographics:","(Baseline=Mostly white)","Non-White: High PSM","Non-White: Low PSM",
          "Overtime work:","(Baseline=Never required)","Occasionally and Frequently required: High PSM","Occasionally and Frequently required: Low PSM",
          "Key job task:","(Baseline=Analysis identifying community needs)","Interaction with community and organizational members: High PSM","Interaction with community and organizational members: Low PSM"))

Full.table <- merge(x=Base.Factor, y=Interaction.Full,
                    by='group', all.x=TRUE)

Full.table$group <- factor(Full.table$group,
                           levels=c("Total pay:", "(Baseline=Slightly below aveage)","About average and Slightly above average: High PSM","About average and Slightly above average: Low PSM",
                                    "Performance bonuses binary:", "(Baseline=No (Fixed salary))","Yes (bonuses): High Self-Efficacy","Yes (bonuses): Low Self-Efficacy",
                                    "Goal based evaluation:","(Baseline=No (Supervisor based)","Yes (goal based): High Self-Efficacy","Yes (goal based): Low Self-Efficacy",
                                    "Current community involvement:","(Baseline=Rare participation)","Moderate and Frequent participation: High PSM","Moderate and Frequent participation: Low PSM",
                                    "Community income:","(Baseline=Low income)","Average and High income: High PSM","Average and High income: Low PSM",
                                    "Community demographics:","(Baseline=Mostly white)","Non-White: High PSM","Non-White: Low PSM",
                                    "Overtime work:","(Baseline=Never required)","Occasionally and Frequently required: High PSM","Occasionally and Frequently required: Low PSM",
                                    "Key job task:","(Baseline=Analysis identifying community needs)","Interaction with community and organizational members: High PSM","Interaction with community and organizational members: Low PSM"))

ggplot(Full.table, aes(x = AME, y = group)) + 
  geom_vline(aes(xintercept = 0), size = .25, linetype = "dashed") + 
  geom_errorbarh(data = Full.table, aes(xmax = upper, xmin = lower), size = .3, height = .6, position = position_nudge(y = 0), color="black") +
  geom_errorbarh(data = Full.table, aes(xmax = upper_90, xmin = lower_90), size = .99, height = .6, position = position_nudge(y = 0), color="black") +
  geom_point(data = Full.table, shape=19, size = 3, position = position_nudge(y = 0), color='black') +
  scale_y_discrete(limits=rev) +
  theme(axis.text.y = element_text(hjust=0)) +
  theme(axis.text.y = element_text(face = 
                                     c("plain","plain","plain","bold","plain","plain","plain","bold","plain","plain","plain","bold",
                                       "plain","plain","plain","bold","plain","plain","plain","bold","plain","plain","plain","bold","plain","plain","plain","bold",
                                       "plain","plain","plain","bold"))) +
  ylab("") +
  xlab("AMCE: Effect on Pr(Government Job Preferred)") +
  theme(axis.text.y = element_text(size=12, colour = "black")) +
  theme(text = element_text(size=12, colour = "black"))


# 2.2.2. Full Model: Continuous measures (Appendix Table F1)
Full.AMCE.Conti <- lm(Chosen_Job ~ `Total pay binary`*PSM_Factor + `Performance bonuses binary`*Efficacy_Factor + `Goal based evaluation`*Efficacy_Factor +
                        `Current community involvement binary`*PSM_Factor + `Community income binary`*PSM_Factor + `Community demographics binary`*PSM_Factor +
                        `Overtime work binary`*PSM_Factor + `Key job task binary`*PSM_Factor, 
                      data=P4P_Clean_Binary)
summary(Full.AMCE.Conti)
Full.AMCE.Conti.cl <- cluster.vcov(Full.AMCE.Conti, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(Full.AMCE.Conti, Full.AMCE.Conti.cl)
Twoway.se.Conti <- sqrt(diag(Full.AMCE.Conti.cl))


Full.AMCE.Collapse.Conti <- lm(Chosen_Job ~ `Total pay`*PSM_Factor + `Performance bonuses`*Efficacy_Factor + `Job performance evaluation`*Efficacy_Factor +
                                 `Current community involvement`*PSM_Factor + `Community income`*PSM_Factor + `Community demographics`*PSM_Factor +
                                 `Overtime work`*PSM_Factor + `Key job task`*PSM_Factor, 
                               data=P4P_Clean_Binary)
summary(Full.AMCE.Collapse.Conti)
Full.AMCE.Collapse.Conti.cl <- cluster.vcov(Full.AMCE.Collapse.Conti, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(Full.AMCE.Collapse.Conti, Full.AMCE.Collapse.Conti.cl)
Twoway.se.Collapse.Conti <- sqrt(diag(Full.AMCE.Collapse.Conti.cl))

stargazer(Full.AMCE.Conti,Full.AMCE.Collapse.Conti,
          se=list(Twoway.se.Conti,Twoway.se.Collapse.Conti), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          order = c(1,2,3,5,6,7,8,10,24,25,26,11,27,28,12,29,30,13,31,32,33,14,34,35,15,36,37,38,
                    4,9,16,39,40,17,41,42,43,18,44,45,46,19,47,48,20,49,50,21,51,52,53,22,54,55,23,56,57,58),
          covariate.labels=c("Total pay binary: About average and Slightly above average",
                             "Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes",
                             "Performance bonuses: large (20%)",
                             "Performance bonuses: moderate (10%)",
                             "Performance bonuses: small (5%)",
                             "Goal based evaluation: Yes",
                             "Goal based evaluation: Attendance",
                             "Goal based evaluation: Changes",
                             "Goal based evaluation: Satisfaction",
                             "Current community involvement binary: Moderate and Frequent participation",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income binary: Average and High income",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics binary: Non-White",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work binary: Occasionally and Frequently required",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task binary: Interaction with community and organizational members",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "PSM (A)",
                             "Self-Efficacy (B)",
                             "Total pay binary: About average and Slightly above average X (A)",
                             "Total pay: About average X (A)","Total pay: Slightly above average X (A)",
                             "Performance bonuses binary: Yes X (B)",
                             "Performance bonuses: large (20%)  X (B)",
                             "Performance bonuses: moderate (10%) X (B)",
                             "Performance bonuses: small (5%) X (B)",
                             "Goal based evaluation: Yes X (B)",
                             "Goal based evaluation: Attendance X (B)",
                             "Goal based evaluation: Changes X (B)",
                             "Goal based evaluation: Satisfaction X (B)",
                             "Current community involvement binary: Moderate and Frequent participation X (A)",
                             "Current community involvement: Frequent participation X (A)",
                             "Current community involvement: Moderate participation X (A)",
                             "Community income binary: Average and High income X (A)",
                             "Community income: Average income X (A)", "Community income: High income X (A)",
                             "Community demographics binary: Non-White X (A)",
                             "Community demographics: Mostly African American X (A)",
                             "Community demographics: Mostly Hispanic X (A)",
                             "Community demographics: Multiracial X (A)",
                             "Overtime work binary: Occasionally and Frequently required X (A)",
                             "Overtime work: Frequently required X (A)","Overtime work: Occasionally required X (A)",
                             "Key job task binary: Interaction with community and organizational members X (A)",
                             "Key job task: Coordination with community groups X (A)",
                             "Key job task: Direct interaction with community residents X (A)",
                             "Key job task: Teamwork with peers and supervisors X (A)"),
          out="Appendix_Table_F1.txt") 



# 2.2.3 PSM Only
# Joint
PSM.AMCE.Joint <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                 `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                 `Overtime work binary` + `Key job task binary` + PSM_C + `Total pay binary`:PSM_C + 
                 `Current community involvement binary`:PSM_C + `Community income binary`:PSM_C + `Community demographics binary`:PSM_C +
                 `Overtime work binary`:PSM_C + `Key job task binary`:PSM_C, 
               data=P4P_Clean_Binary)
summary(PSM.AMCE.Joint)
PSM.AMCE.Joint.cl <- cluster.vcov(PSM.AMCE.Joint, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(PSM.AMCE.Joint, PSM.AMCE.Joint.cl)
Twoway.se.Joint <- sqrt(diag(PSM.AMCE.Joint.cl))

# Binary
PSM.AMCE.Totalpay <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                       `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                       `Overtime work binary` + `Key job task binary` + PSM_C + `Total pay binary`:PSM_C, 
                     data=P4P_Clean_Binary)
summary(PSM.AMCE.Totalpay)
PSM.AMCE.Totalpay.cl <- cluster.vcov(PSM.AMCE.Totalpay, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(PSM.AMCE.Totalpay, PSM.AMCE.Totalpay.cl)
Twoway.se.Totalpay <- sqrt(diag(PSM.AMCE.Totalpay.cl))

PSM.AMCE.Cominv <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                          `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                          `Overtime work binary` + `Key job task binary` + PSM_C + `Current community involvement binary`:PSM_C, 
                        data=P4P_Clean_Binary)
summary(PSM.AMCE.Cominv)
PSM.AMCE.Cominv.cl <- cluster.vcov(PSM.AMCE.Cominv, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(PSM.AMCE.Cominv, PSM.AMCE.Cominv.cl)
Twoway.se.Cominv <- sqrt(diag(PSM.AMCE.Cominv.cl))

PSM.AMCE.Cominc <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                        `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                        `Overtime work binary` + `Key job task binary` + PSM_C + `Community income binary`:PSM_C, 
                      data=P4P_Clean_Binary)
summary(PSM.AMCE.Cominc)
PSM.AMCE.Cominc.cl <- cluster.vcov(PSM.AMCE.Cominc, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(PSM.AMCE.Cominc, PSM.AMCE.Cominc.cl)
Twoway.se.Cominc <- sqrt(diag(PSM.AMCE.Cominc.cl))

PSM.AMCE.Comdem <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                        `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                        `Overtime work binary` + `Key job task binary` + PSM_C + `Community demographics binary`:PSM_C, 
                      data=P4P_Clean_Binary)
summary(PSM.AMCE.Comdem)
PSM.AMCE.Comdem.cl <- cluster.vcov(PSM.AMCE.Comdem, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(PSM.AMCE.Comdem, PSM.AMCE.Comdem.cl)
Twoway.se.Comdem <- sqrt(diag(PSM.AMCE.Comdem.cl))

PSM.AMCE.Overtime <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                        `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                        `Overtime work binary` + `Key job task binary` + PSM_C + `Overtime work binary`:PSM_C, 
                      data=P4P_Clean_Binary)
summary(PSM.AMCE.Overtime)
PSM.AMCE.Overtime.cl <- cluster.vcov(PSM.AMCE.Overtime, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(PSM.AMCE.Overtime, PSM.AMCE.Overtime.cl)
Twoway.se.Overtime <- sqrt(diag(PSM.AMCE.Overtime.cl))

PSM.AMCE.Jobtask <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                          `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                          `Overtime work binary` + `Key job task binary` + PSM_C + `Key job task binary`:PSM_C, 
                        data=P4P_Clean_Binary)
summary(PSM.AMCE.Jobtask)
PSM.AMCE.Jobtask.cl <- cluster.vcov(PSM.AMCE.Jobtask, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(PSM.AMCE.Jobtask, PSM.AMCE.Jobtask.cl)
Twoway.se.Jobtask <- sqrt(diag(PSM.AMCE.Jobtask.cl))

# Reporting (Appendix Table F2)
A <- PSM.AMCE.Joint
B <- PSM.AMCE.Totalpay
C <- PSM.AMCE.Cominv
D <- PSM.AMCE.Cominc
E <- PSM.AMCE.Comdem
F <- PSM.AMCE.Overtime
G <- PSM.AMCE.Jobtask
stargazer(A,B,C,D,E,F,G,
          se=list(Twoway.se.Joint,Twoway.se.Totalpay,Twoway.se.Cominv,Twoway.se.Cominc,
                  Twoway.se.Comdem,Twoway.se.Overtime,Twoway.se.Jobtask), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay binary: About average and Slightly above average",
                             "Performance bonuses binary: Yes",
                             "Goal based evaluation: Yes",
                             "Current community involvement: Moderate and Frequent participation",
                             "Community income: Average and High income",
                             "Community demographics: Non-White",
                             "Overtime work: Occasionally and Frequently required",
                             "Key job task: Interaction with community and organizational members",
                             "High PSM (A)",
                             "Total pay binary: About average and Slightly above average X (A)",
                             "Current community involvement: Moderate and Frequent participation X (A)",
                             "Community income: Average and High income X (A)", 
                             "Community demographics: Non-White X (A)",
                             "Overtime work: Occasionally and Frequently required X (A)",
                             "Key job task: Interaction with community and organizational members X (A)"),
          out="Appendix_Table_F2.txt") 


# 2.2.3 Self-Efficacy
# Joint
SE.AMCE.Joint <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                       `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                       `Overtime work binary` + `Key job task binary` + Efficacy_C + `Performance bonuses binary`:Efficacy_C +
                      `Goal based evaluation`:Efficacy_C, 
                     data=P4P_Clean_Binary)
summary(SE.AMCE.Joint)
SE.AMCE.Joint.cl <- cluster.vcov(SE.AMCE.Joint, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(SE.AMCE.Joint, SE.AMCE.Joint.cl)
Twoway.se.Joint <- sqrt(diag(SE.AMCE.Joint.cl))

# Binary
SE.AMCE.Bonus <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                      `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                      `Overtime work binary` + `Key job task binary` + Efficacy_C + `Performance bonuses binary`:Efficacy_C, 
                    data=P4P_Clean_Binary)
summary(SE.AMCE.Bonus)
SE.AMCE.Bonus.cl <- cluster.vcov(SE.AMCE.Bonus, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(SE.AMCE.Bonus, SE.AMCE.Bonus.cl)
Twoway.se.Bonus <- sqrt(diag(SE.AMCE.Bonus.cl))

SE.AMCE.Goal <- lm(Chosen_Job ~ `Total pay binary` + `Performance bonuses binary` + `Goal based evaluation` +
                      `Current community involvement binary` + `Community income binary` + `Community demographics binary` +
                      `Overtime work binary` + `Key job task binary` + Efficacy_C + `Goal based evaluation`:Efficacy_C, 
                    data=P4P_Clean_Binary)
summary(SE.AMCE.Goal)
SE.AMCE.Goal.cl <- cluster.vcov(SE.AMCE.Goal, cbind(P4P_Clean_Binary$ID, P4P_Clean_Binary$Treatment.ID))
coeftest(SE.AMCE.Goal, SE.AMCE.Goal.cl)
Twoway.se.Goal <- sqrt(diag(SE.AMCE.Goal.cl))

# Reporting (Appendix Table F3)
stargazer(SE.AMCE.Joint,SE.AMCE.Bonus,SE.AMCE.Goal,
          se=list(Twoway.se.Joint,Twoway.se.Bonus,Twoway.se.Goal), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay binary: About average and Slightly above average",
                             "Performance bonuses binary: Yes",
                             "Goal based evaluation: Yes",
                             "Current community involvement: Moderate and Frequent participation",
                             "Community income: Average and High income",
                             "Community demographics: Non-White",
                             "Overtime work: Occasionally and Frequently required",
                             "Key job task: Interaction with community and organizational members",
                             "High Self-Efficacy (A)",
                             "Performance bonuses binary: Yes X (A)",
                             "Goal based evaluation: Yes X (A)"),
          out="Appendix_Table_F3.txt") 



# 2.3 Non-linearity
P4P_Interflex <- P4P_Clean_Binary %>%
  rename(Totalpaybinary=`Total pay binary`,
         Performancebonusesbinary=`Performance bonuses binary`,
         Goalbasedevaluation=`Goal based evaluation`,
         Currentcommunityinvolvementbinary=`Current community involvement binary`,
         Communityincomebinary=`Community income binary`,
         Communitydemographicsbinary=`Community demographics binary`,
         Overtimeworkbinary=`Overtime work binary`,
         Keyjobtaskbinary=`Key job task binary`)

#install.packages('interflex', type = "source", repos = 'http://cran.us.r-project.org')
library(interflex)

# Binning estimator
TotalPay <- interflex(Y = "Chosen_Job", D = "Totalpaybinary", X = "PSM_Factor",
                 Z = c("Performancebonusesbinary","Goalbasedevaluation","Currentcommunityinvolvementbinary",
                       "Communityincomebinary","Communitydemographicsbinary","Overtimeworkbinary",
                       "Keyjobtaskbinary"),
                 data = P4P_Interflex, estimator = "binning", 
                 base = "Slightly below average", cl=c("ID","Treatment.ID"),
                 theme.bw = TRUE)
TotalPay.f <- plot(TotalPay,ylab="Total pay",xlab="PSM")
print(TotalPay$tests$p.wald)

Bonus <- interflex(Y = "Chosen_Job", D = "Performancebonusesbinary", X = "Efficacy_Factor",
                 Z = c("Totalpaybinary","Goalbasedevaluation","Currentcommunityinvolvementbinary",
                       "Communityincomebinary","Communitydemographicsbinary","Overtimeworkbinary",
                       "Keyjobtaskbinary"),
                 data = P4P_Interflex, estimator = "binning", 
                 base = "No (fixed salary)", cl=c("ID","Treatment.ID"),
                 theme.bw = TRUE)
Bonus.f <- plot(Bonus,ylab="Performance bonuses",xlab="Self-Efficacy")
print(Bonus$tests$p.wald)

Goal <- interflex(Y = "Chosen_Job", D = "Goalbasedevaluation", X = "Efficacy_Factor",
                   Z = c("Totalpaybinary","Performancebonusesbinary","Currentcommunityinvolvementbinary",
                         "Communityincomebinary","Communitydemographicsbinary","Overtimeworkbinary",
                         "Keyjobtaskbinary"),
                   data = P4P_Interflex, estimator = "binning", 
                   base = "No (supervisor based)", cl=c("ID","Treatment.ID"),
                   theme.bw = TRUE)
Goal.f <- plot(Goal,ylab="Goal based evaluation",xlab="Self-Efficacy")
print(Goal$tests$p.wald)

Cominv <- interflex(Y = "Chosen_Job", D = "Currentcommunityinvolvementbinary", X = "PSM_Factor",
                      Z = c("Totalpaybinary","Performancebonusesbinary","Goalbasedevaluation",
                            "Communityincomebinary","Communitydemographicsbinary","Overtimeworkbinary",
                            "Keyjobtaskbinary"),
                      data = P4P_Interflex, estimator = "binning", 
                      base = "Rare participation", cl=c("ID","Treatment.ID"),
                      theme.bw = TRUE)
Cominv.f <- plot(Cominv,ylab="Community involvement",xlab="PSM")
print(Cominv$tests$p.wald)

Cominc <- interflex(Y = "Chosen_Job", D = "Communityincomebinary", X = "PSM_Factor",
                    Z = c("Totalpaybinary","Performancebonusesbinary","Goalbasedevaluation",
                          "Currentcommunityinvolvementbinary","Communitydemographicsbinary","Overtimeworkbinary",
                          "Keyjobtaskbinary"),
                    data = P4P_Interflex, estimator = "binning", 
                    base = "Low income", cl=c("ID","Treatment.ID"),
                    theme.bw = TRUE)
Cominc.f <- plot(Cominc,ylab="Community income",xlab="PSM")
print(Cominc$tests$p.wald)

Comdem <- interflex(Y = "Chosen_Job", D = "Communitydemographicsbinary", X = "PSM_Factor",
                    Z = c("Totalpaybinary","Performancebonusesbinary","Goalbasedevaluation",
                          "Currentcommunityinvolvementbinary","Communityincomebinary","Overtimeworkbinary",
                          "Keyjobtaskbinary"),
                    data = P4P_Interflex, estimator = "binning", 
                    base = "Mostly white", cl=c("ID","Treatment.ID"),
                    theme.bw = TRUE)
Comdem.f <- plot(Comdem,ylab="Community demographics",xlab="PSM")
print(Comdem$tests$p.wald)

Overtime <- interflex(Y = "Chosen_Job", D = "Overtimeworkbinary", X = "PSM_Factor",
                    Z = c("Totalpaybinary","Performancebonusesbinary","Goalbasedevaluation",
                          "Currentcommunityinvolvementbinary","Communityincomebinary","Communitydemographicsbinary",
                          "Keyjobtaskbinary"),
                    data = P4P_Interflex, estimator = "binning", 
                    base = "Never required", cl=c("ID","Treatment.ID"),
                    theme.bw = TRUE)
Overtime.f <- plot(Overtime,ylab="Overtime work",xlab="PSM")
print(Overtime$tests$p.wald)

Task <- interflex(Y = "Chosen_Job", D = "Keyjobtaskbinary", X = "PSM_Factor",
                      Z = c("Totalpaybinary","Performancebonusesbinary","Goalbasedevaluation",
                            "Currentcommunityinvolvementbinary","Communityincomebinary","Communitydemographicsbinary",
                            "Overtimeworkbinary"),
                      data = P4P_Interflex, estimator = "binning", 
                      base = "Analysis identifying community needs", cl=c("ID","Treatment.ID"),
                      theme.bw = TRUE)
Task.f <- plot(Task,ylab="Key Job Task",xlab="PSM")
print(Task$tests$p.wald)

# Plot (Appendix Figure F2)
library(cowplot)
P4P_Mod <- cowplot::plot_grid(TotalPay.f,Bonus.f,Goal.f,Cominv.f,
                              Cominc.f,Comdem.f,Overtime.f,Task.f,
                              ncol=2, nrow=4)
P4P_Mod


# Kernel estimator
TotalPay.k <- interflex(Y = "Chosen_Job", D = "Totalpaybinary", X = "PSM_Factor",
                      Z = c("Performancebonusesbinary","Goalbasedevaluation","Currentcommunityinvolvementbinary",
                            "Communityincomebinary","Communitydemographicsbinary","Overtimeworkbinary",
                            "Keyjobtaskbinary"),
                      data = P4P_Interflex, estimator = "kernel", 
                      base = "Slightly below average", cl=c("ID","Treatment.ID"),
                      theme.bw = TRUE)
TotalPay.k.f <- plot(TotalPay.k,ylab="Total pay",xlab="PSM")

Bonus.k <- interflex(Y = "Chosen_Job", D = "Performancebonusesbinary", X = "Efficacy_Factor",
                   Z = c("Totalpaybinary","Goalbasedevaluation","Currentcommunityinvolvementbinary",
                         "Communityincomebinary","Communitydemographicsbinary","Overtimeworkbinary",
                         "Keyjobtaskbinary"),
                   data = P4P_Interflex, estimator = "kernel", 
                   base = "No (fixed salary)", cl=c("ID","Treatment.ID"),
                   theme.bw = TRUE)
Bonus.k.f <- plot(Bonus.k,ylab="Performance bonuses",xlab="Self-Efficacy")

Goal.k <- interflex(Y = "Chosen_Job", D = "Goalbasedevaluation", X = "Efficacy_Factor",
                  Z = c("Totalpaybinary","Performancebonusesbinary","Currentcommunityinvolvementbinary",
                        "Communityincomebinary","Communitydemographicsbinary","Overtimeworkbinary",
                        "Keyjobtaskbinary"),
                  data = P4P_Interflex, estimator = "kernel", 
                  base = "No (supervisor based)", cl=c("ID","Treatment.ID"),
                  theme.bw = TRUE)
Goal.k.f <- plot(Goal.k,ylab="Goal based evaluation",xlab="Self-Efficacy")

Cominv.k <- interflex(Y = "Chosen_Job", D = "Currentcommunityinvolvementbinary", X = "PSM_Factor",
                    Z = c("Totalpaybinary","Performancebonusesbinary","Goalbasedevaluation",
                          "Communityincomebinary","Communitydemographicsbinary","Overtimeworkbinary",
                          "Keyjobtaskbinary"),
                    data = P4P_Interflex, estimator = "kernel", 
                    base = "Rare participation", cl=c("ID","Treatment.ID"),
                    theme.bw = TRUE)
Cominv.k.f <- plot(Cominv.k,ylab="Community involvement",xlab="PSM")

Cominc.k <- interflex(Y = "Chosen_Job", D = "Communityincomebinary", X = "PSM_Factor",
                    Z = c("Totalpaybinary","Performancebonusesbinary","Goalbasedevaluation",
                          "Currentcommunityinvolvementbinary","Communitydemographicsbinary","Overtimeworkbinary",
                          "Keyjobtaskbinary"),
                    data = P4P_Interflex, estimator = "kernel", 
                    base = "Low income", cl=c("ID","Treatment.ID"),
                    theme.bw = TRUE)
Cominc.k.f <- plot(Cominc.k,ylab="Community income",xlab="PSM")

Comdem.k <- interflex(Y = "Chosen_Job", D = "Communitydemographicsbinary", X = "PSM_Factor",
                    Z = c("Totalpaybinary","Performancebonusesbinary","Goalbasedevaluation",
                          "Currentcommunityinvolvementbinary","Communityincomebinary","Overtimeworkbinary",
                          "Keyjobtaskbinary"),
                    data = P4P_Interflex, estimator = "kernel", 
                    base = "Mostly white", cl=c("ID","Treatment.ID"),
                    theme.bw = TRUE)
Comdem.k.f <- plot(Comdem.k,ylab="Community demographics",xlab="PSM")

Overtime.k <- interflex(Y = "Chosen_Job", D = "Overtimeworkbinary", X = "PSM_Factor",
                      Z = c("Totalpaybinary","Performancebonusesbinary","Goalbasedevaluation",
                            "Currentcommunityinvolvementbinary","Communityincomebinary","Communitydemographicsbinary",
                            "Keyjobtaskbinary"),
                      data = P4P_Interflex, estimator = "kernel", 
                      base = "Never required", cl=c("ID","Treatment.ID"),
                      theme.bw = TRUE)
Overtime.k.f <- plot(Overtime.k,ylab="Overtime work",xlab="PSM")

Task.k <- interflex(Y = "Chosen_Job", D = "Keyjobtaskbinary", X = "PSM_Factor",
                  Z = c("Totalpaybinary","Performancebonusesbinary","Goalbasedevaluation",
                        "Currentcommunityinvolvementbinary","Communityincomebinary","Communitydemographicsbinary",
                        "Overtimeworkbinary"),
                  data = P4P_Interflex, estimator = "kernel", 
                  base = "Analysis identifying community needs", cl=c("ID","Treatment.ID"),
                  theme.bw = TRUE)
Task.k.f <- plot(Task.k,ylab="Key Job Task",xlab="PSM")

# Plot (Appendix Figure F3)
P4P_Mod.k <- cowplot::plot_grid(TotalPay.k.f,Bonus.k.f,Goal.k.f,Cominv.k.f,
                              Cominc.k.f,Comdem.k.f,Overtime.k.f,Task.k.f,
                              ncol=2, nrow=4)
P4P_Mod.k



# 2.4 Separate sample analysis (Three-way Interactions)
table(P4P_Clean_Binary$Study.ID,P4P_Clean_Binary$Race_B,exclude=NULL)

P4P_Clean_Binary_Race <- P4P_Clean_Binary %>% subset(Race_B!="Prefer not to answer")

Full.AMCE.Two <- lm(Chosen_Job ~ `Total pay binary`*Race_B + `Performance bonuses binary`*Race_B + `Goal based evaluation`*Race_B +
                        `Current community involvement binary`*Race_B + `Community income binary`*Race_B + `Community demographics binary`*Race_B +
                        `Overtime work binary`*Race_B + `Key job task binary`*Race_B, 
                      data=P4P_Clean_Binary_Race)
summary(Full.AMCE.Two)
Full.AMCE.Two.cl <- cluster.vcov(Full.AMCE.Two, cbind(P4P_Clean_Binary_Race$ID, P4P_Clean_Binary_Race$Treatment.ID))
Full.stat.Two <- coeftest(Full.AMCE.Two, Full.AMCE.Two.cl)
Two.se <- sqrt(diag(Full.AMCE.Two.cl))

Full.AMCE.Three <- lm(Chosen_Job ~ `Total pay binary`*PSM_C*Race_B + `Performance bonuses binary`*Efficacy_C*Race_B + `Goal based evaluation`*Efficacy_C*Race_B +
                        `Current community involvement binary`*PSM_C*Race_B + `Community income binary`*PSM_C*Race_B + `Community demographics binary`*PSM_C*Race_B +
                        `Overtime work binary`*PSM_C*Race_B + `Key job task binary`*PSM_C*Race_B, 
                      data=P4P_Clean_Binary_Race)
summary(Full.AMCE.Three)
Full.AMCE.Three.cl <- cluster.vcov(Full.AMCE.Three, cbind(P4P_Clean_Binary_Race$ID, P4P_Clean_Binary_Race$Treatment.ID))
Full.stat.Three <- coeftest(Full.AMCE.Three, Full.AMCE.Three.cl)
Three.se <- sqrt(diag(Full.AMCE.Three.cl))

# Reporting (Appendix Table F4)
stargazer(Full.AMCE.Two,Full.AMCE.Three,
          se=list(Two.se,Three.se), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          order = c(1,4,6,7,8,9,10,11,
                    2,5,3,12,15,18,20,22,24,26,28,
                    13,16,19,21,23,25,27,29,
                    14,17,30,31,32,33,34,35,36,37),
          covariate.labels=c("Total pay binary: About average and Slightly above average",
                             "Performance bonuses binary: Yes",
                             "Goal based evaluation: Yes",
                             "Current community involvement binary: Moderate and Frequent participation",
                             "Community income binary: Average and High income",
                             "Community demographics binary: Non-White",
                             "Overtime work binary: Occasionally and Frequently required",
                             "Key job task binary: Interaction with community and organizational members",
                             "PSM (A)",
                             "Self-Efficacy (B)",
                             "Minority (C)",
                             "Total pay binary: About average and Slightly above average X (A)",
                             "Performance bonuses binary: Yes X (B)",
                             "Goal based evaluation: Yes X (B)",
                             "Current community involvement binary: Moderate and Frequent participation X (A)",
                             "Community income binary: Average and High income X (A)",
                             "Community demographics binary: Non-White X (A)",
                             "Overtime work binary: Occasionally and Frequently required X (A)",
                             "Key job task binary: Interaction with community and organizational members X (A)",
                             "Total pay binary: About average and Slightly above average X (C)",
                             "Performance bonuses binary: Yes X (C)",
                             "Goal based evaluation: Yes X (C)",
                             "Current community involvement binary: Moderate and Frequent participation X (C)",
                             "Community income binary: Average and High income X (C)",
                             "Community demographics binary: Non-White X (C)",
                             "Overtime work binary: Occasionally and Frequently required X (C)",
                             "Key job task binary: Interaction with community and organizational members X (C)",
                             "PSM (A) X Minority (C)",
                             "Self-Efficacy (B) X Minority (C)",
                             "Total pay binary: About average and Slightly above average X (A) X (C)",
                             "Performance bonuses binary: Yes X (B) X (C)",
                             "Goal based evaluation: Yes X (B) X (C)",
                             "Current community involvement binary: Moderate and Frequent participation X (A) X (C)",
                             "Community income binary: Average and High income X (A) X (C)",
                             "Community demographics binary: Non-White X (A) X (C)",
                             "Overtime work binary: Occasionally and Frequently required X (A) X (C)",
                             "Key job task binary: Interaction with community and organizational members X (A) X (C)"),
          out="Appendix_Table_F4.txt")

